Decay of scalar variance in isotropic turbulence in a bounded domain 
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Abstract. The decay of scalar variance in isotropic turbulence in a bounded domain is investi- 
gated. Extending the study of Touil, Bertoglio and Shao [J. Turbul. 03:49 2002] to the case of a 
passive scalar, the effect of the finite size of the domain on the lengthscales of turbulent eddies and 
scalar structures is studied by truncating the infrared range of the wavenumber spectra. Analytical 
arguments based on a simple model for the spectral distributions show that the decay exponent 
for the variance of scalar fluctuations is proportional to the ratio of the Kolmogorov constant to 
the Corrsin-Obukhov constant. This result is verified by closure calculations in which the Corrsin- 
Obukhov constant is artificially varied. Large-Eddy Simulations provide support to the results and 
give an estimation of the value of the decay exponent and of the scalar to velocity time scale ratio. 



I. INTRODUCTION 



Most turbulent flows on earth are bounded by walls. 
An important part of theoretical studies of turbulence is 
however devoted to homogeneous fields, ruling out the ex- 
istence of boundaries that are present in most situations 
of practical interest or laboratory experiments. Even in 
nearly homogeneous flows, such as decaying grid turbu- 
lence, the integral lengthscale will grow indefinitely so 
that, if the turbulence was initially strong enough, this 
lengthscale will sooner or later become comparable to 
the size of the domain. Evidently, the understanding of 
spatially bounded turbulence is of primary importance 
in a wide range of academic and engineering problems 
and remains one of the major challenges in turbulence 
research. 

Indeed, the presence of walls adds an enormous com- 
plexity to the description of turbulence. Without taking 
into account the zoology of structures and phenomena 
associated with the existence of a strong shear near the 
boundaries present in most wall bounded flows, there is a 
simple academic case where global predictions can easily 
be performed. This is the case of an isotropic turbulence 
freely evolving in a bounded domain. Following Tennekes 
and Lumley [l] it can be predicted that, after a period of 
free decay, in which the kinetic energy evolves classically, 
as soon as the integral lengthscale becomes limited by 
the size of the domain, the turbulent kinetic energy k de- 
cays as t~ 2 (and therefore that the dissipation e and rms 
vorticity decay as t~ 3 and as t~ 3 ^ 2 respectively). This 
prediction was confirmed experimentally by Skrbek and 
Stalp Q in confined superfluid grid turbulence. After the 
lengthscale saturation occured, they measured a t~ 3 / 2 
decay for the rms vorticity and performed an analysis of 
the results based on a simple model spectrum with an 
infrared cut-off. It was already anticipated by Bertoglio 
and Jeandel 0] , that the scale limitation in spectral stud- 
ies of turbulence could be roughly taken into account by 
imposing an infrared cut-off in the energy spectrum, rep- 
resenting the fact that eddies larger than the domain size 
cannot exist. Direct numerical simulations (DNS), Large- 



Eddy Simulations (LES) and two-point closure calcula- 
tions [4] confirmed the result for isotropic turbulence. 
The effect of a bounded domain on anisotropic turbu- 
lence was similarly investigated by Biferale et aL[5j by 
DNS. The present work continues this line of research 
by investigating the decay of passive scalar fluctuations 
in a confined turbulent flow. This effort contributes to 
a better understanding of turbulent mixing in bounded 
domains. 

As stated above, the decay exponent of the kinetic en- 
ergy can be estimated by using the relation (Tennekes 
and Lumley [l[) 
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with d the size of the domain that determines the integral 
lengthscale. Assuming self-similar decay of the turbulent 
energy spectrum, power law decay can be expected for k 
and e: 

k~{t-t )- n e~n(;-< )- ( " +1) (2) 

with to a virtual origin. Using ((T|) and the relation 

k t = -e, (3) 

one obtains 

(t-t )-^ ~dr 1 (t-t )-i n (4) 

so that the equality of the exponents leads to the already 
mentioned result n — 2. 

The equation for the variance of passive scalar fluc- 
tuations 2 in isotropic turbulence without mean scalar 
gradients is 
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,2 ,t = -e e . 



(5) 



in which eg is the dissipation of scalar fluctuations eg. 

It is tempting to apply the above reasoning for the 
kinetic energy to the decay of the passive scalar. The 
equivalent of |T]) for scalar decay is: 
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FIG. 1: Model spectrum used to predict selfsimilar decay of 
turbulence. 



assuming power laws for k, eg and 6 2 yields: 

(t - t r {ne+1) ~ d-\t - io )-(™«+«/ 2 ) (7) 

which yields by substituting n — 2 and using the equality 
of the exponents 



(n e + 1) = -{rig + 1) 



(8) 



which is verified for every ng. We therefore need to use 
another method to determine ng. This is the purpose 
of section [ill m which an analytical study based on a 
simplified spectral form of the energy and scalar spectra 
is performed that allows to derive analytical expressions 
for the scalar decay. These expressions are then validated 
by two-point closure computations of the eddy-damped 
quasi-normal Markovian (EDQNM) type and confirmed 
by Large-Eddy Simulations in section IIIII All through 
the paper the Schmidt number will be assumed unity. 



II. ANALYSIS AND SCALINGS 

For the purpose of clarity some results for the velocity 
field are recalled, before starting the discussion of the pas- 
sive scalar. The decay of kinetic energy in an unbounded 
isotropic turbulent field can be roughly predicted by con- 
sidering the simple model spectrum (c.f. Comte-Bcllot 
and Corrsin Q), 

AK S for K < K L 

E(K) = I C K e 2 / 3 K-^ 3 for K L < K < K n (9) 







for K> K„ 



where Ck is the Kolmogorov constant. This spectrum 
corresponds to the sketch in figure [T] It will be assumed 
in the following that K n ^> Kl- E(K) is related to the 
kinetic energy k by 



k = 



E(K)dK 



(10) 



Assuming power laws for k and e and that A is con- 
stant during the time evolution, the following decay ex- 
ponent is obtained 

2(s + l) , s 

n = v J . (11) 



Note that the constancy of A is a reasonable assumption 
for s < 4 (see Lesieur and Schertzer Q). A study of the 
time evolution of A for s = 4 can be found for example in 
Chasnov 0], but is outside the scope of the present paper. 
The present work focuses on the behavior of turbulence 
once the growth of the integral scales is limited by the 
domain size. As shown in Skrbek and Stalp 2] or Touil 
et al. [J] , as soon as the integral scale becomes comparable 
to the domain size, the decay exponent tends to 2, a 
value larger than before saturation (n = 10/7 for s = 4 
for example). This exponent 2 can be deduced assuming 
an energy spectrum with a sharp infrared cut-off, of the 
form 

'0 for K < K mf 

E{K) = { C K e 2 / 3 K- 5 / 3 for K inf < K < K v (12) 







for K > Kn 



i.e. an inertial range extending from the wavenumber 
Ki n f = 27r/d, with d the domain size, to K n . 

Integrating (|12|) gives in the limit of an infinite 
Reynolds number: 



(13) 



Assuming power law decay for e ((2|) and using ([3]) one 
can express e as: 



C\ (n + 1) 3 

■ (t - t )- (n+1) 



(14) 



so that the equality between time exponents in the bal- 
ance implies that n = 2. 

Let us now come back to the central issue of this work, 
the decay of the variance of passive scalar fluctuations. 
By analogy with ^ the scalar spectrum can be assumed 
to have the form: 

' A K S ' for K< Kg 

Eg(K) = { C CO e- 1/3 egK- 5 / 3 for Kg < K < K n 







for K > K n 



(15) 

with Ceo the Corrsin-Obukhov constant. The scalar 
variance is related to the spectrum by 



i— r°° 

-9 2 = / Eg(K)dK. 
z Jo 



(16) 



Let us assume that 9 2 decays as a power law: 

¥^(t-t y na (17) 

then eg decays as 

e e ~n 9 (t-t )-( n ° +1 l (18) 

This is not a trivial assumption and it will be checked 
a posteriori by closure computations. The decay of 9 2 



3 



depends not only on its spectral distribution, but also on 
the injection scale of the scalar fluctuations 0, [kJ Ell- 
in the present work we consider the case in which the 
initial scales of the scalar field and the velocity field are 
comparable, Kl ~ Kg. Within this framework, power 
law decay can be possible fl2l ]. and an expression simi- 
lar to pip can be derived for the scalar variance. The 
resulting expression (Herring et aZ. (l3j) is: 



ng 



2{s' + l) 



(19) 



An important quantity directly related to the decay of 
scalar variance and appearing in most engineering mod- 
els, is the velocity to scalar time scale ratio, defined by: 



2ke g 



(20) 



The assumption of power law decay of the integral quan- 
tities ((2]) and (|17| yields immediately the result r = ng/n, 
which gives for the present case 



s'+l 

8 + 1' 



(21) 



A more elaborate expression, depending on the Schmidt 
number, was recently proposed by Ristorcelli [llj |. 

To investigate the effect of the limitation of the tur- 
bulent and scalar lengthscales when the integral length- 
scales become comparable to the size of the domain, we 
still assume that the energy spectrum is given by (112|) 
and we furthermore postulate that the scalar spectrum 
has the form 

'0 for K < K inf 

Eg{K) = { C CO e- 1/3 eeK- 5 / 3 for K lnf < K < K n 
for K > K n . 

(22) 

Using expression (| 14[) for the dissipation in (|22|) and in- 
tegrating to obtain 9 2 : 

l^r 3 



d2— r> a ^-i/s^- 2 / 3 



Hnf > ( 2 3) 

and using (fT7j) together with ((5]) gives 



" no%¥-(t - t r ( " e+1) = -(* ~ t )- {ns+1) (24) 

ZCk 

so that equality of the exponents, leads again to expres- 
sion © which is satisfied whatever the value of ng . To 
deduce an estimation of ng we must therefore use the 
equality between the prefactors in (f2"4"|) . It is found that 



ng = 2 



Ceo 



The time scale ratio directly follows 
r = ng/2 



(25) 



(26) 



It must be emphasized that the situation is different from 
what is found for the kinetic energy spectrum, or for the 
scalar freely decaying in an unbounded domain, where 
the values of the decay exponents are found by simply 
using the relation between the exponents. The fact that 
here use is made of a relation between the prefactors, 
leads to the rather untypical situation where the Corrsin- 
Obukhov and Kolmogorov constants expicitly appear in 
the expression of the decay exponents. This is for exam- 
ple in contrast with what is found for the velocity field 
where n is not a function of the Kolmogorov constant, 
but is entirely determined by the power law exponent s 
of the spectrum in the low wavenumber region. 

The relations (|25|) and (|26|) stress the importance of an 
precise knowledge of Ck and Ceo- A review aimed at 
determining the precise value of Ceo by a compilation of 
numerous experimental results was performed by Sreeni- 
vasan [3]. In this work the value of C' co , the constant 
intervening in the one-dimensional scalar spectrum was 
found to be approximately 0.3 < C' co < 0.6. C\ 



J co 



is 



related to Cqo by the relation Ceo — (5/3) C' co which 
gives 0.5 < Ceo < 1-0 according to the experiments. 

Using EDQNM and LES it will be investigated in the 
following section if the predictions of r and ng are right. 



III. CLOSURE AND LARGE-EDDY 
SIMULATION 

A. EDQNM 

In this paper computations are performed within the 
framework of the eddy-damped quasi-normal Marko- 
vian (EDQNM) closure, as proposed_for isotropic tur- 
bulence by Orszag 15] and Leith [l6|. The isotropic 
scalar formulation used here was fir st p roposed by Vi- 
gnon et a/.[l7], HH and Herring et a/. [HI]. The evolution 
equations for the kinetic energy spectrum of an isotropic 
turbulence and for the passive scalar spectrum are re- 
spectively: 



-2vK 2 E{K) + T{K) 



d_ 

Ot 



Eg(K) = ~2nK 2 Eg(K) + Tg(K) 



(27) 



(28) 



Time dependence is omitted for notational simplicity. 
The viscosity v is equal to the diffusivity k for the case 
of unity Schmidt number, which we consider here. The 
non-linear transfer terms T{K) and Tg(K) are expressed 
using the classical EDQNM formulation for isotropic tur- 
bulence and scalar field: 



t{k) = f e KPQ {t) 

Ja(k) 



xy + z c 

Q 



E(Q) x 



K 2 E(P) - P 2 E{K)^dPdQ (29) 
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FIG. 2: Time evolution of the kinetic energy spectrum (top) 
and scalar variance spectrum (bottom) 



T 9 (K)= f 6 6 K pQ (t)^-^-E(Q) x 
Ja(k) 



Q 3 



K 3 E 9 {P) - KP 2 Eg(K) \dPdQ (30) 



in which the damping coeficients are expressed using the 
classical forms: 



r)(K) = xJ J R 2 E(R)dR 



R 2 E(R)dR 



n"(K) = \ 2 J R 2 E(R)dR 



(33) 



Three constants have to be fixed. For A we use the clas- 
sical value A = 0.355. The values for Ai and A2 re- 
quire more attention. We follow Herring et aZ. [l3| and 
choose a zero value for Ai for compatability with the 
Lagrangian History Direct Interaction Approximation of 
Kraichnan The value of A2 is not fixed and is re- 

lated to the Corrsin-Obukhov constant. We will first use 
the classical value (see Lesieur [l2| and Herring et al.fu§) 
A 2 = 1.3. The results will be analyzed in the following 
section. Then, we will take advantage of the fact that 
varying A2 provides a simple and efficient way to vary 
the Corrsin-Obukhov constant in the model and we will 
use this degree of freedom to check relations (|2"5j) and 

The computational domain ranges from JQ n / (low- 
wavenumber or infrared cut-off, related to the size d 
of the bounded domain by i"Q„/ = 27r/d) to 4K V , K v 
being the Kolmogorov wavenumber). The resolution is 
approximately 14 wavenumbers per decade. The initial 
spectrum is a von Karman spectrum (20| . The energy- 
containing range is characterized by wavenumber Kl , the 
wavenumber at which E(K) has its maximum. At time 
t = 0, initial conditions are such that Kl is greater than 
Ki n f. The Taylor-scale Reynolds number evaluated at 
the time when Kr, ~ Ki n j is approximately 10 5 . This 
high value for the Reynolds number is chosen to allow 
a precise determination of the Kolmogorov and Corrsin- 
Obukhov constants. 



in which A denotes the domain such that the three wave 
vectors , P, Q form a triangle and x, y and z are the 
cosines of the angles respectively opposite to K, P, Q in 
this triangle. The EDQNM characteristic times are given 
by 



OkpqH) 



I _ e -(u(K 2 + P 2 +Q 2 )+ n {K)+ri(P)+v{Q))t 

u(K 2 + P 2 + Q 2 ) + j](K) + r)(P) + n(Q) 

(31) 



B. Large- Eddy Simulation 

The code used for the LES computations is a classi- 
cal pseudo-spectral code with fourth-order Runge-Kutta 
time integration scheme. The resolution is 128 3 grid 
points. As initial spectrum, a von Karman spectrum 
is used which behaves as K 4 at small K and as K~ 5 ^ 3 
for large K. The CZZS dynamic model [HI [13 is used 
to model the subgrid stress and scalar flux. 



IV. RESULTS 



7 KPQ 



(*) 



1 _ e -(MK 2 + P 2 )+"Q 2 W(K)+ri'(P)+ri"(Q))t 

k(K 2 + P 2 ) + vQ 2 + rf(K) + n'(P) + if'(Q) 

(32) 



In Figure [5] EDQNM results are shown for the time 
evolution of the kinetic energy spectrum and the scalar 
variance spectrum. It can be observed that both spectra 



5 



display a clear K 5 / 3 inertial range. The inertial ranges 
are observed before saturation and are still present after. 
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FIG. 3: Decay of kinetic energy and scalar fluctuations. 
EDQNM computations results 




FIG. 4: Compensated kinetic energy spectrum (left) and scalar 
variance spectrum (right) 



These results were obtained with the classical values 
of the constants (A = 0.355, Ai = and A2 = 1.3). The 
kinetic energy and scalar variance are shown in figure [3] 
It is observed that after a short time the kinetic energy 
decays following a power law with an exponent close to 



the classical value 10/7. As soon as the large scales be- 
come saturated by the presence of the infrared cut-off, at 
t sa t, the power law exponent changes and takes a value 
of 2 as predicted in the preceding sections. At very long 
times, the behavior changes, corresponding to a final pe- 
riod of viscous decay. It is observed in figure [3j that the 
scalar variance also decays following a power law with an 
exponent close to 10/7. At longer times, after saturation 
has occurred, the decay is found to follow a i~ 3 - 25 power 
law. The Kolmogorov and Corrsin-Obukhov constants 
are determined by carefully evaluating the compensated 
spectra. Figure 2] shows the compensated spectra E(K) 
and Eg(K). The Kolmogorov constant is Ck ~ 1-35, 
a slightly low value. The Corrsin-Obukhov constant is 
then approximately 0.85. The ratio 2Ck/Cco yields 3.2 
which is close to the value of ne observed in figure [3] 

Figure [5] shows the time evolution of the velocity to 
scalar time scale ratio. It is observed that after a tran- 
sient, r goes to a constant value close to 1 corresponding 
to the freely decaying period where both k and 6 2 evolve 
as t~ w / 7 . At longer times, after saturation has occurred, 
a second plateau is observed with a value close to 1.6. 
This value is in agreement with the analysis in section HT1 
that predicts that r = Ck/Cco- 




FIG. 5: Time evolution of the velocity to scalar time scale 
ratio 



In figure [6] the kinetic energy and scalar variance from 
the Large-Eddy Simulation are shown. For the saturated 
decay a power law exponent for the energy close to 2 is 
observed in agreement with the results of Touil et al. [1] . 
The scalar variance decays with an exponent close to 4. 
This corresponds according to (|2"5")l to a ratio Ck/Cco 
close to 2. This is well verified as can be oberved in 
figure [7] where the ratio Ck/Cco is plotted as a function 
of K/K c with K c the cut-off wavenumber. 

To check the relations (|23|) and (|26[) in more detail 
Ceo is now artificially varied by changing the value of 
A2 in the range [0.7-2.0]. By examining the compen- 
sated spectra it was observed that this corresponded to a 
Corrsin-Obukhov constant that takes values in the range 
[0.46-1.28]. The time evolutions of 8 2 are shown in fig- 
ure [5] for three different values of A2. Clear power laws 
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FIG. 6: Decay of kinetic energy and scalar fluctuations. LES 
results 
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FIG. 7: Ratio Ck/Cco from the Large-Eddy Simulations 



results of the closure and the analytical results (|25|) and 
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FIG. 9: Scalar decay exponent ng as a function of the ratio of 
the Kolmogorov Ck and Corrsin-Obukhov constant Ceo (top). 
Scalar decay exponent as a function of the time scale ratio r 
(bottom) 



are observed. The power law exponent is affected by the 
variation of A2 . 
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FIG. 8: Decay of passive scalar variance for varying A2 



In figure M we test the results (I25t and against 
our computations. Good agreement is found between the 



V. CONCLUSION 

Analytical expressions were derived for the decay ex- 
ponent of the scalar variance and the velocity to scalar 
time scale ratio for a passive scalar decaying in bounded 
isotropic turbulence. These expressions depend on the 
Kolmogorov and Corrsin-Obukhov constants. It is rather 
uncommon in turbulence theory that a decay exponent 
depends on the values of these inertial range constants. 
The proposed relations were tested against EDQNM 
computations in which the Corrsin-Obukhov constant 
was artificially varied. Good agreement was observed 
and clear power law decay was observed for time evolu- 
tion of the scalar variance. The values for the decay ex- 
ponent and Kolmogorov and Corrsin-Obukhov constants 
obtained by Large-Eddy Simulation are consistent with 
these results. The decay exponent is shown to be approx- 
imately 4 in the LES, a significantly larger value than the 
exponent for the turbulent kinetic energy which in this 
case is close to 2. The present study shows the impor- 
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tance of a precise knowledge of spectral constants as their 
values can directly determine the temporal behavior of 
integral quantities. 

An experimental verification of the results presented 
in this paper would probably be difficult to achieve, as 
in usual laboratory experiments, the Reynolds number 
is too low for the wall bounded regime to be observed 
before the final viscous decay occurs. To overcome this 
difficulty, Skrbek and Stalp [2J performed measurements 
in Helium superfluid. They were able to measure the 
—3/2 decay exponent for the rms vorticity that corre- 



sponds to the bounded domain regime, but no data on 
scalar fluctuations were reported in their experimental 
study. 
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